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In this paper we analyze the second-order Kuramoto model presenting a positive correlation be¬ 
tween the heterogeneity of the connections and the natural frequencies in scale-free networks. We 
numerically show that discontinuous transitions emerge not just in disassortative but also in assor¬ 
tative networks, in contrast with the first-order model. We also find that the effect of assortativity 
on network synchronization can be compensated by adjusting the phase damping. Our results show 
that it is possible to control collective behavior of damped Kuramoto oscillators by tuning the 
network structure or by adjusting the dissipation related to the phases movement. 
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I. INTRODUCTION 

Synchronization is pervasive in nature, society and 
technology [1]. This collective behavior emerges from 
the interaction of neurons in the central nervous system, 
power grids, crickets, heart cells and lasers [mg. Syn- 
chronization arises due to the adjustment of rhythms of 
self-sustained periodic oscillators weakly connected CHS 
and can be treated mathematically by the model pro¬ 
posed by Kuramoto [3j. The general Kuramoto model 
assumes that the natural frequencies of the oscillators 
are selected from unimodal and symmetric random dis¬ 
tributions [T]. In this case, a second-order phase transi¬ 
tion to synchronization can be observed mm- However, 
the first-order Kuramoto model can exhibit discontinu¬ 
ous phase transitions HE]- For instance, in one of the 
first works on this topic, Pazo [4] showed that, if uniform 
frequency distributions are considered, first-order tran¬ 
sitions emerge in fully connected Kuramoto oscillators. 
Very recently, Gomez-Gardehes et al. [7) verified that a 
discontinuous synchronization transition can also occur 
in scale-free networks as an effect of positive correlation 
between the natural frequencies and network topology. 
This discovery has triggered many ensuing works, which 
analyzed explosive synchronization analytically and nu¬ 
merically [81416]. 

The early works on explosive synchronization (e.g. & 
doi) suggested that the correlation between frequency 
and degree distributions is the only condition required for 
the emergence of a discontinuous synchronization tran¬ 
sition in scale-free networks. However, subsequent pa- 
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pers have shown that different criteria to set the fre¬ 
quency mismatch between the oscillators m, the pres¬ 
ence of time-delay [18], non-vanishing degree-degree cor¬ 
relation dadllEg or the inclusion of noise m can dra¬ 
matically change the type of the phase transitions, even 
in the regime of fully connected graphs [13 ng. 

In the case of degree-degree correlation, Li et al. ca 
verified that assortative scale-free networks no longer un¬ 
dergo a discontinuous transition, even if the network 
presents a positive correlation between structural and dy¬ 
namical properties. This behavior was also observed in 
the synchronization of FitzHugh-Nagumo (FHN) oscilla¬ 
tors coupled in scale-free networks under the constraint 
of correlating frequencies and degrees eg. Furthermore, 
Zhu et al. m found that discontinuous transitions only 
emerge in networks subjected simultaneously to nega¬ 
tive degree-degree and frequency-frequency correlations. 
Effects of degree-degree correlation on general network 
synchronization phenomena were also analyzed in litera¬ 
ture (cf. P3JE5] for studies in the context of the master 
stability function formalism). For instance, Bernardo et 
al. [21| studied scale-free networks of identical Rossler os¬ 
cillators and showed that disassortative mixing enhances 
network synchronization, when compared with uncorre¬ 
lated networks m On the other hand, regarding the 
synchronization of weighted networks, assortative mixing 
can enhance synchronization, depending on the weight¬ 
ing procedure [23] . 

All these works considered only the first-order Ku¬ 
ramoto model. However, in various systems, the second- 
order Kuramoto model is more suitable to describe the 
emergence of synchronization [^177] . Indeed, many real- 
world systems, such as power-grid networks [28JI35] . su¬ 
perconducting Josephson junctions [26| and many other 
applications [211123, can be modeled as networks whose 
units are second-order Kuramoto oscillators. 

In the context of explosive synchronization, the second- 
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order model was introduced only recently [Mi 157] . In 
contrast to [7], where the authors verified that nodes in 
scale-free network join the synchronous group abruptly 
giving rise to a discontinuous synchronization transition, 
we have shown that this behavior is no longer observed 
when an inertia term is included ISSUST]. Despite also 
observing a discontinuous transition of the order parame¬ 
ter, we verified that in the second-order Kuramoto model 
with frequencies proportional to degrees, nodes join the 
synchronous component successively grouped into clus¬ 
ter of nodes with the same degree, a phenomenon called 
cluster explosive synchronization (CES) |36j . 


The influence of network structure on the emergence of 
explosive synchronization in the second-order Kuramoto 
model proposed in [36] has not been addressed yet, since 
only uncorrelated networks have been considered [MUSH- 
Among important network properties, the degree-degree 
correlation is observed in several complex networks |38h 
HD] and it plays a fundamental role in many dynamical 
processes, such as epidemic spreading and synchroniza¬ 
tion ED For example, the degree-degree correlation can 
change the type of phase transitions of the first-order 
Kuramoto model with a positive correlation between fre¬ 
quency and degree distributions [m US]. In this way, 
since the frequency mismatch between oscillators has 
been shown to play a crucial role in the emergence of 
abrupt transitions in the first-order Kuramoto model, it 
is natural to ask about the effects of degree-degree cor¬ 
relations on the overall dynamics in models with inertia. 


In this paper we study the second-order Kuramoto 
model in networks with degree-degree correlations, i.e., 
non-vanishing assortativity. We find that the synchro¬ 
nization diagrams have a strong dependence on the net¬ 
work assortativity, but in a different fashion compared to 
the first-order model m- In fact, for negative and posi¬ 
tive assortativity values, the synchronization is observed 
to be discontinuous, depending on the damping coeffi¬ 
cient. Moreover, the upper branch in the synchronization 
diagrams associated to the case in which the coupling 
is decreased is barely affected by different assortativity 
values, again in contrast with the first-order Kuramoto 
model [HHni nzi [la EOj- In other words, we show here 
that one is able to control the hysteretic behaviour of 
the second-order Kuramoto model by tuning the net¬ 
work properties, the phenomenon that was not investi¬ 
gated before. In order to compare with different choices 
of frequencies distributions, we also investigate the dy¬ 
namics of damped Kuramoto oscillators in assortative 
networks using unimodal and even distributions, with¬ 
out being correlated with the local topology. Similarly 
as in the case of frequencies proportional to degrees, we 
again observe very similar behavior for the onset of syn¬ 
chronization over networks with different degree-degree 
correlations. 


II. SYNCHRONIZATION IN CORRELATED 
NETWORKS 

We consider networks where each node is a phase os¬ 
cillator evolving according to the second-order Kuramoto 
model [Il|3| 

^2 f) ^ 

~di = + sin m) - eM , (1) 

i=i 

where a is the dissipating parameter, A is the coupling 
strength and fti is the natural frequency of oscillator 
i, defined according to a given probability distribution 
The heterogeneity of the network connections is 
accounted by the adjacency matrix A = {Aij}, whose 
elements Aij = 1 if oscillators i and j are connected, and 
Aij = 0 otherwise. The collective dynamics of the oscil¬ 
lators is measured by the macroscopic order parameter, 
defined as 

1 ^ 

=-Y^e^OAt) ^ ( 2 ) 

i=i 

where the modulus 0 < r(t) < 1 and '^(t) is the av¬ 
erage phase of the oscillators. The system governed by 
Eqs. exhibits hysteretic synchrony [23112]. The onset 
of synchronization (r > 0) is characterized by a critical 
coupling AJ. when the coupling strength is progressively 
increased from a given Aq. On the other hand, starting 
at synchronouos state and decreasing progressively the 
coupling strength, the oscillators fall into an incoherent 
state (r ^ 0) at coupling A^ < Aj, [27ll4^ . 

Here we study the second-order Kuramoto model (see 
Eq. in which the natural frequency distribution ^(O) 
is correlated with the degree distribution P{k) as Vti = 
ki — (k) [36], where ki is the degree of the oscillator i and 
(k) is the average degree of the network. At first glance 
that particular choice for the frequency assignment could 
sound odd, however it is not difficult to find physical 
scenarios where this configuration is plausible. For ex¬ 
ample, such correlation between dynamics and network 
topology can arise as a consequence of a limited amount 
of resources or energy supply for the oscillators. In fact, 
studies on optimization of synchronization in complex 
networks [431145] have shown that, for a given fixed set 
of allowed frequencies {f2i, ^42, •••, the configuration 
that maximizes the network synchronization is reached 
for cases in which frequencies are positively correlated 
with degrees. Therefore, this correlation between fre¬ 
quencies and local topology can be seen as an optimal 
scenario for the emergence of collective behavior in com¬ 
plex networks. 

We study networks presenting non-vanishing degree- 
degree correlation. Such a correlation is quantified by 
a measure known as assortativity coefficient, A , which 
is the Pearson coefficient between degrees at the end of 
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FIG. 1: Synchronization diagram r(A) with (a) a = 0.2, 
(b) OL — 0.4, (c) OL — 0.6 and (d) ol — 0.8 for assortativity 
A — —0.3. With increasing a, onset of synchronization and 
hysteresis decreases. The natural frequency of each oscillator 
is Q^i — ki — {k) and the networks have N = 10^ and (k) =6. 
The degree distribution follows a power-law P{k) ~ k~^, 
where 7 = 3. Curves in which points are connected by solid 
lines (resp. dashed lines) correspond to the forward (resp. 
backward) continuations of the coupling strength A. 


FIG. 2: Synchronization diagram r(A) for (a) A = —0.3, (b) 
A = —0.1, (c) A = 0.1 and (d) A = 0.3. The dissipation 
coefficient is fixed at a = 1 with natural frequencies given 
by Q^i — ki — (k), as in Fig. All networks considered have 
N = 10^, (k) = 6 and P{k) ~ k~^, where 7 = 3. 


III. NUMERICAL RESULTS 


each link [38], i.e. 

jih - [M-i E, 1 7 + h)] ^ 

Ei\{fi + kl) - [M-i Ei Uji + ki )]' ’ 

where —l<A<l,ji and ki are the degrees associated 
to the two ends of the edge i {i = 1,...,M) and M is 
the total number of edges in the network. In order to 
tune the degree of assortativity of each network, we use 
the method proposed in [46]. The algorithm allows us 
to obtain networks with a desired value of assortativity 
without changing the degree of each node. At each step, 
two edges are selected at random and the four nodes as¬ 
sociated to these edges are ordered from the lowest to 
the highest degree. In order to produce assortative mix¬ 
ing (Al > 0), with a probability p, one new edge connects 
the first and the second node and another new edge links 
the third and fourth nodes. In the case when one of the 
two new edges already exists, the step is discarded and 
a new pair of edges are chosen. This same heuristic can 
also generate disassortative networks (Al < 0) with only 
a slight change in the algorithm. After selecting the four 
nodes and sorting them with respect to their degrees, 
one must rewire, with probability p, the highest degree 
node with the lowest one and, likewise, the second and 
third nodes. After rewiring the network, if the degree of 
assortativity is higher or smaller than the designed Al, 
p is decreased or increased respectively and the network 
is rewired following the procedures described above. In 
order to avoid dead loops, the increasing and decreasing 
steps of p should not be equally spaced. 


In this section we present the results obtained by nu¬ 
merically evolving the equations of motions considering 
Eq. in assortative networks constructed according to 
the model described in the previous section. In all simu¬ 
lations the initial networks are constructed through the 
Barabasi-Albert (BA) model with (k) = 6 and N = 
1 X 10^. 

The order parameter r is calculated with forward and 
backward continuations of the coupling strength A. More 
specifically, by increasing the value of A adiabatically, 
we integrate the system long enough and calculate the 
stationary value of r for each coupling Aq, Aq + ^A, ..., Ao + 
nSX. Similarly, for the backward continuation, we start 
at the value X = n6X^ Xq and decrease A by amounts of 
6X until A = Aq. In both processes we use 6X = 0.5. 

We investigate the dependency of the hysteresis on the 
dissipation parameter a. Fig. shows the forward and 
backward synchronization diagrams r(A) for networks 
with assortativity A = —0.3, but different values of a 
within the interval [0.2,1]. As we can see, the area of 
hysteresis and the critical coupling for the onset of syn¬ 
chronization in the increasing branch tends to decrease 
as a is increased, which also contributes to increase the 
maximal value of the order parameter. 

Next, we fix the dissipation coefficient a = 1 and 
vary the network assortativity in the interval [—0.3,0.3]. 
Fig. i shows the synchronization diagram r(A) for net¬ 
works with different values of assortativity. As A in¬ 
creases, the hysteresis becomes less clear and the on¬ 
set of synchronization in the decreasing branch tends 
to increase. Surprisingly, the critical coupling of the in- 
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FIG. 3: (Color online) Contour plot on a-A plane colored according to (a) the maximal order parameter MO, (b) the maximal 
order parameter difference MV and (c) the hysteresis area S. A and a are varied within the interval [—0.3, 0.3] and [0.2,1], 
respectively. For each pair {a, A), 40 times simulations are performed with the coupling strength in the interval as in Fig. H 
i.e, A G [0, 25]. 


creasing branch for the second-order Kuramoto model is 
weakly affected, which is in sharp contrast with results 
concerning models without inertia [13 EOl- More pre¬ 
cisely, in the first-order Kuramoto model with frequen¬ 
cies correlated with degrees, the critical coupling for the 
onset of synchronization in scale-free networks increases 
as the network becomes more assortative Hi Ho]. The 
same phenomenon was observed in the synchronization 
of FHN oscillators [12]. 

In order to evaluate more accurately the dependency 
of the synchronization transitions on the level of assor- 
tativity and the dissipation parameter, we introduce the 
maximal order parameter MO, the maximal order pa¬ 
rameter difference MV and the hysteresis area S in the 
synchronization diagrams as a function of A and a, re¬ 
spectively, as follows, 

MO = (max (r(A))), (4) 



FIG. 4: Parameters are the same as in Fig.[^ except that the 
natural frequencies are randomly selected from a Lorentzian 
distribution g(V) = + 11^)) . 


A4P = (max |r^(A) — r^(A)l) , (5) 

5=(y \T\\)-T^{\)\d\y (6) 

where A G [Aq, Aq + nSX], (.) denotes the average of dif¬ 
ferent realizations and | . | the absolute value. r^(A) and 
r^(A) are the order parameters for increasing and de¬ 
creasing coupling strength A, respectively. If r(A) in¬ 
creases as A grows, then MO is usually obtained at 
the maximal coupling strength Aq + nSX, i.e. MO = 
r(Ao +n(5A). MV G [0,1] quantifies the hysteresis differ¬ 
ence. If the system shows a continuous phase transition 
with a perfect match between increasing and decreasing 
coupling strength diagrams, then 5^0. In this way, S is 
a quantitative index to evaluate the hysteretic behavior 
related to the emergence of synchronization. 

Comparing Figs.andwe observe a clear dependence 
of the synchronization diagrams on the assortativity A 
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FIG. 5: Parameters are the same as in Fig. except that 
natural frequencies are randomly selected from a Lorentzian 
distribution 
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FIG. 6: (Color online) Contour plots similar as in Fig. for networks with a natural frequency distribution given by g{Q) = 
l/(7r(l + n2)). 


and on the dissipation parameter a. Moreover, note that 
it is also possible to obtain similar dependencies of r on 
A by selecting different values of a and A. In order to 
better grasp this apparent equivalence in the dynamical 
behavior of the system for different choices of the param¬ 
eters A and a, we show in Fig. [^the quantities defined 
in Eqs. [^and as a function of a and A. As we can 
see in Fig.Wa), similar values for the maximal order pa¬ 
rameter AiO are obtained according to the initial setup 
of the model. More specifically, the level of synchroniza¬ 
tion of the network can be chosen by tuning the assor- 
tativity or the dissipation parameter in the dynamical 
model. Therefore, for the second-order Kuramoto model 
in the case of frequencies positively correlated with de¬ 
gree, high levels of coherent behavior are obtained by ei¬ 
ther strongly assortative or disassortative networks, once 
the dissipation parameter a is properly selected. Inter¬ 
estingly, the maximal gap between the increasing and 
decreasing branches quantified by MV (Fig. [^b)) has a 
maximum around A — —0.1 and a G [0.4,0.^, showing 
that the area of metastability in the stability diagram of 
the model [31] is maximized for this set of parameters. 

A similar effect can also be observed for the hystere¬ 
sis area S in Fig. [^c). The maximal S in the synchro¬ 
nization diagram is reached for networks with A = —0.1 
and a = 0.2. Furthermore, similar values of S are ob¬ 
tained by different sets of a and M, which shows an inter¬ 
esting interplay between the topological parameter (as- 
sortativity) and the dynamical one (dissipation). More 
precisely, topological properties related to degree-degree 
correlations can be counterbalanced by the dissipation 
parameter in the dynamical model. This property could 
have interesting applications in the control of synchro¬ 
nization in networks modeled by the second-order Ku¬ 
ramoto model. In particular, if one is interested to reduce 
hysteresis in a system, such task can be accomplished by 
either increasing the dissipation or the degree mixing in 
the network. Therefore, the question usually addressed 
in studies regarding the first-order Kuramoto model that 
is whether assortativity could enhance synchronization 
or not [mill] turns out to be harder to answer for the 
damped version of the model. The reason for that is that 
the asymptotic behavior of the system strongly depends 


on the combination of parameters A and a, which al¬ 
lows at the same time much more options to control the 
system by tuning such parameters. 

In order to analyze how assortative mixing influences 
the dynamics of networks of damped Kuramoto oscil¬ 
lators without the constraint of having Qi oc ki^ we 
also compute the same forward and backward synchro¬ 
nization diagrams considering a Lorentzian distribution 
g{V) = for different values of degree assortativ¬ 

ity A. Similarly as before, as a first experiment, we fix 
the assortativity at A = —0.3 and vary the dissipation 
parameter a as indicated in Fig.|^ Again, as we increase 
a the hysteresis area tends to decrease. The same effect 
is observed for a fixed a with varying M, as depicted in 
Fig.|5] 

Calculating MO as a function of a and A we note, 
however, a slight different dependence compared to the 
case where frequencies are proportional to degree. As we 
can see in Fig.j^a), for a fixed value of a, MO is weakly 
affected by the change of the degree mixing, except for 
the case a = 0.2. Nonetheless, the model with frequen¬ 
cies correlated with degrees presents larger fluctuations 
for the maximum value of coherent behavior, comparing 
Figs.l^a) andl^a). Furthermore, the maximum value of 
MV ^ig.j^b)) for Lorentzian frequency distributions is 
obtained for slight disassortative networks with low val¬ 
ues of the dissipation parameter a. We also note that, 
for large a, the maximum gap MV starts to decrease, in 
contrast to the case with frequencies correlated with de- 
gree (Fig.|^b)), where MV is close to zero for almost the 
entire range considered of a for which M > 0. Finally, an¬ 
alyzing S in Fig.j^c) we note the same interplay between 
topological perturbations in the networks, accounted by 
changes in assortative mixing, and dynamical features in 
the oscillator model characterized by the dissipation pa¬ 
rameter a. As shown in this figure, similar values of S 
are achieved by controlling the parameters a and M, and 
highly and poorly hysteretic synchronization diagrams 
can be obtained by different strategies, i.e., changing the 
network structure (M) or the dynamical nature of the 
oscillators (a). 
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IV. CONCLUSIONS 

First-order synchronization transitions for the Ku- 
ramoto model in complex networks have been known as 
a consequence of positive correlation between network 
structure, represented by the degree distribution, and the 
intrinsic oscillatory dynamics, represented by the natural 
frequency distribution of the oscillators [7Hin]. 

Here, we have numerically shown that such transitions 
for the second-order Kuramoto model also depend on the 
degree mixing in the network connection. More precisely, 
discontinuous transitions of networks of second-order Ku¬ 
ramoto oscillators can take place not only in disassorta- 
tive ones but also in assortative ones, in contrast to what 
has been observed for the first-order Kuramoto model 
in which the correlation between topology and dynamics 
is also present m- The reason behind this phenomenon 
can be regarded as an effect of the dynamical equivalence 
of changes in the network structure, played by assortative 
mixing, and changes in the oscillator model (dissipation 
parameter). In other words, a given final configuration 
of a network of second-order Kuramoto oscillators can 
be achieved by tuning the network structure or by ad¬ 
justing the dissipation related to the phases movement. 
As previously mentioned, this finding can have impor¬ 
tant applications on controlling network synchronization 
where, for instance, there are costs associated to lead the 
system to a given desirable state, allowing the adoption 


of different strategies to accomplish such task. 

Our results show that the hysteretic behavior of the 
order parameter vanishes for some assortativities, sug¬ 
gesting that the transition might become continuous. 
However, to properly determine the nature of the phase 
transition, a more detailed study should be addressed. 
Moreover, the theoretical description of sychronization in 
correlated networks is still an open problem and mean- 
field theories that account for degree-degree correlations 
should also be developed. As a future study, it would 
be interesting to further analyze the present model re¬ 
lating the recent approaches on mean-field approxima¬ 
tion of first-order Kuramoto oscillators in assortative 
networks m and the low-dimensional behavior of the 
second-order model [48] . 
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